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Abstract. 

Chimera states—curious symmetry-broken states in systems of identical coupled 
oscillators—typically occur only for certain initial conditions. Here we analyze 
their basins of attraction in a simple system comprised of two populations. Using 
perturbative analysis and numerical simulation we evaluate asymptotic states and 
associated destination maps, and demonstrate that basins form a complex twisting 
structure in phase space. Understanding the basins’ precise nature may help in the 
development of control methods to switch between chimera patterns, with possible 
technological and neural system applications. 


PACS numbers: 05.45.-a, 05.45.Xt, 05.65.+b 


Keywords: chimera states, basins of attraction, hierarchical network, neural networks. 


Submitted to: New J. Phys. 



Basins of Attraction for Chimera States 

1. Introduction 


2 


Self-emergent synchronization is a key process in networks of coupled oscillators, and 
is observed in a remarkable range of systems, including pendulum clocks, pedestrians 
on a bridge locking their gait, Josephson junctions, flashing fireflies, the beating of 
the heart, circadian clocks in the brain, chemical oscillations, metabolic oscillations 
in yeast, life cycles of phytoplankton, and genetic oscillators [1-13]. About a decade 
ago, a study [14] revealed the existence of chimera states, in which a population of 
identical coupled oscillators splits up into two parts, one synchronous and the other 
incoherent. This state is counter-intuitive as it appears despite the oscillators being 
identical. Recent experiments using metronomes, (electro-)chemical oscillators and 
lasing systems [15-19] have demonstrated the existence of chimera states in real-world 
settings; previous theoretical studies have also confirmed the robustness of chimeras 
subjected to a range of adverse conditions, including additive noise, varied oscillator 
frequencies, varied coupling topologies, and other imperfections [20-30]. 

Chimeras are known to arise in systems with nonlocal coupling that decays 
with increasing distance between phase oscillators, thus bridging the gap between the 
extremes of local (nearest-neighbor) and global (all-to-all) coupling];. Such long-range 
coupling is characteristic of many real-world technological and biological [34-36] systems. 
In many systems, chimeras are steady-state solutions stably coexisting with the fully 
synchronized stateS, not emerging via spontaneous symmetry breaking, and are thus 
only attained via a certain class of initial conditions [14,24,29]. Figure 1 graphically 
demonstrates this puzzling aspect of basins of attraction for chimera states: apparently 
similar initial conditions (panel B) can evolve to completely different steady-states (panel 
C). Thus, a natural question arising in any practical situation is: given a random initial 
phase configuration, how likely is the system to converge to a chimera state? Even 
though this important question was raised in 2010 [39], basins of attraction for chimera 
states have not yet been investigated systematically. 

One difficulty with the examination of basins of attractions is that they are 
computationally expensive to obtain, e.g. via Monte Carlo simulation [40]. Here, in 
contrast, we use primarily analytic methods to explain the structure of the phase space 
and provide a systematic study of the basins of attraction leading to chimeras in the 
thermodynamic limit. 

Model. The simplest realization of nonlocal coupling is achieved with two 
populations, where each population is more strongly coupled to itself than to the 
neighboring population (see Figure 1 panel A). It has been used as a model for several 
investigations of chimera states [15,26,27,41-44]; here chimeras manifest themselves 

]; For oscillators with non-constant amplitudes it appears that local and global coupling can both be 
sufficient [31-33]. 

S Known exceptions where chimera states and fully synchronized states are not necessarily bi-stable 
are certain generalizations of the present system with amplitude dynamics [33,37] or non-linear delay 
feedback [21,38]. 
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Figure 1. A: Schematic of the system under investigation. B: Three superficially 
similar oscillator phase distributions taken as initial conditions. C: Oscillator phase 
distributions after long-time evolution of system—each corresponds to the initial 
condition shown directly above it. DS: “desync-sync” state; SSq: “sync-sync” state; 
“SD”: “sync-desync” state. 


as a state with one synchronous and one asynchronous population. Accordingly, we 
consider the Kuramoto-Sakaguchi model with n = 2 populations [41,42] each of size 
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ei - a). 


( 1 ) 


cr'=l 1=1 

where Of is the phase of the A:th oscillator k = in population a e {1,2} 

and u! is the oscillator frequency. For consistency with previous work [41,42,45], 
we assume the coupling is symmetric with neighbor-coupling iFo-o-' = A'cr'o- = ^ and 
self-coupling = /i. Imposing without loss of generality fi + n = 1, the coupling 
can be parameterized by the coupling disparity A = — n. We redefine the phase 
lag parameter via (3 = t :/2 — a as chimeras emerge in the limit of near-cosine¬ 
coupling {P —7- 0) for this type of system [41,45,46]. The mean field order parameter 

(iOj) describes the synchronization level of population a 
with Rfj = 1 for perfect and < 1 for partial synchronization. We consider the 
thermodynamic limit N'^ oo, allowing us to express the ensemble dynamics in terms of 
the continuous oscillator density /‘^(d,C(j). This facilitates a low-dimensional description 
of the dynamics via the Ott-Antonsen (OA) ansatz [47-49] in terms of the mean-field 
order parameter of each population, = f f‘^{0,t)d0 with 0 < Pa < 1, 

see Appendix A and Appendix B. 

By virtue of the translational symmetry const., the resulting dynamics 
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Figure 2. State variables (pi, P 2 ) are interpreted as cylindrical coordinates. Phase 
space is structured by (i) two invariant rays, Rq and (dashed); and (ii) two invariant 
surfaces, Si and S 2 , forming the side and top surfaces of the cylinder slab. Except for a 
set of measure zero, all trajectories converge to one of three locations: SD chimera state 
on Si (red), DS chimera state on S 2 (blue), or fully synchronized state SSq (yellow). 
(A,/3) = (0.1,0.025); filled/empty circles denote stable/unstable fixed points. Small 
yellow dots denote initial conditions. 


are effectively three dimensional with the angular phase difference — 4 > 2 ^ obeying 

i-p? 


pi = 
P2 = 
Ip = 


2 

l-pj 

2 

1 + pj 

2p2 

1 + pf 

2pi 


[l^pi sin /5 + i'p 2 sin (/3 — 

[pp 2 sin /5 + upi sin (/3 + 

[pp 2 cos /3 + upi cos (/5 + 

[ppi cos P + up 2 cos {P — ip )], 


( 2 ) 

(3) 

(4) 


on domain D = {(pi, P 2 , ip)\0 < pi ,2 < 1, —ti" <ip< ti"}- 

Phase space is visualized using cylindrical coordinates (pi, P 2 ), see Figure 2. The 
translation II : /3 i->- /3 + tt reverses time in Eqs. (2)-(4), thus inverting flow in phase 
space and stability of fixed points; we restrict our attention to 0 < /? < tt in what 
follows. 


2. Invariant Manifolds and fixed points 

Analysis of equations (2)-(4) reveals the existence of two invariant surfaces defined by 
So- = {(pi, P 2 , t^)|P(T = 1} C n (the top (blue) and lateral (red) surfaces of the cylinder 





















Basins of Attraction for Chimera States 


5 


displayed in Figure 2). The dynamics on these manifolds were studied previously [41]: 
chimera states are born in a saddle-node bifurcation and undergo a Hopf bifurcation for 
larger coupling disparity A = /i — u. The resulting stable limit cycle grows with A until 
eventually it is destroyed in a homoclinic bifurcation. Studying the basins of attraction, 
we generalize the previous analysis by considering the entire three-dimensional phase 
space D. 

Numerically, we observe that all trajectories with initial conditions pi,p 2 < 1 
are attracted to one of the invariant surfaces. From there, any of three attractors 
can be asymptotically approached: (i) a partially synchronized limit point (stable 
chimera, either SD or DS), (ii) a limit cycle (breathing chimera, either SD or DS), 
or (hi) the fully synchronized state SSq at {pi,P 2 ,'f’) = (1,1,0). Furthermore, unstable 
fixed points exist: a fully synchronized state SStt at {pi,p 2 ,'ip) = (1,1, vr), and several 
unstable saddle chimeras (UC) (see [28] and Appendix D). The dynamics on Si and 
S 2 are related due to the invariance of Eqs. (2)-(4) under the symmetry operation 
S : (pi,p2,'0) ^ {p2,pi,-'ip). 

Outside of Si and S 2 , trajectories follow a complex winding motion, structured 
around the two invariant rays Rq and R.,^ defined by pi = p 2 with ^ = 0 and ip = tt, 
respectively (see Figure 2 and Appendix C). Other than the origin, which is a repeller, 
there are no fixed points for pi, p 2 < 1 (see Appendix D). Thus, limit cycles in the interior 
of the phase space are also absent. In principle, a chaotic attractor could appear inside 
D but is not observed. 

3. Numerical investigations 

First insights regarding basins of attraction for chimera states were gathered via 
simple Monte Carlo integration of uniformly distributed random initial conditions for 
Pi,P 2 G [0,1] and Ip G [—7r,7r]. These computations reveal that the probability p{A, (3) 
of ending up in a chimera state depends primarily on (3 with a maximum value for 
/5 —>■ 0, see Figure 3. This approach provides information about the sizes of the basins 
of attraction, but it reveals little about their structure. We therefore ask: how is the 
three-dimensional phase space structured? 

To better reflect symmetries of the phase space, the dynamics may be re-expressed 
in terms of the sum and difference of the order parameters (see Figure 2), s = \{pi+ P 2 ) 
with 0 < s < 1, and d= \{pi — P 2 ), with —a(s) < d < a{s) where a{s) = | — s| (see 

Eqns. (B.6)-(B.8)). 

In order to characterize the structure of the basins of attraction, we compute the 
destination maps for a set of initial conditions (s, d, ip). Eigure 4A shows a typical cross- 
section of the destination map with fixed s: basins form a spiraling structure around 
(the ray {d,ip) = (0,7r)), with SD and DS basins always separated by the (often 
thin) basin for SSq. The thickness of the basin spiral arms increases away from R^^, with 
maximum near Rq (the ray {d,ip) = (0,0)). 

The area ratio between basins for SD (or, by symmetry, DS) and SSq is related 
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Figure 3. Probabilities to obtain chimera states via random sampling of initial 
conditions {pi,p 2 ,ip). Chimeras appear within the wedge defined by a saddle-node 
bifurcation (SN, solid) for small A and a homoclinic bifurcation (HC, dotted) for large 
A [41]. Phase portraits in the p 2 = I plane are shown (insets) with stable nodes 
(full circles), unstable chimera (white circle) and saddle chimera (half-filled circles), 
together with its stable (solid) and unstable (dashed) manifolds. For intermediate 
A, the asynchronous order parameter undergoes a Hopf bifurcation (HB, dashed). 
Probabilities for ending up in either SD/DS chimera were measured by realizing 1000 
random initial conditions {pi,p 2 ,ip) for each parameter value set. Further details are 
in Appendix E. 


to the probability that a random initial condition will lead to a chimera state, and 
depends on parameters A and P as follows. For P ^ 0, the SSq basin occupies an 
infinitesimal fraction of the area. As P increases, the SSq basin increases its area until it 
occupies the entire plane at P = /5sn(A) when the chimera state is annihilated through 
a saddle-node bifurcation. For A < Asn(/ 5) or A > Ahc(/ 5), no chimera state exists and 
the entire basin belongs to the SSq state. With increasing A > Asn, the (total) basin 
area of SSq gradually decreases from 100% approaching a constant near the homoclinic 
bifurcation, see Figure E2. 

As s increases from zero, basins merge and pinch-off in an alternating fashion 
(see Figure 4, Sec. 4 and Supplementary Video 1) so that the basin boundaries rotate 
clockwise about Rq {{d,'tp) = (0,0) in Figure 4A). Once s reaches Sc ~ Vl — A, this 
rotation stops, demonstrating that knowledge of the trajectory position in the s = Sc 
plane is sufficient for determining its final fate. 

The basin density appears singular near R.,^, with a nested structure that allows 
even tiny perturbations of the initial condition to strongly influence the final state. More 
generally, the highly alternating basin structure is reflected in the times to reach steady- 
state attractors, which are displayed in Figure 4 B. Figures 5 and El show destination 
times T along a section of that figure between the origin and {d,'tp) = (0,7r), revealing 
a power-law behavior that may result from this nested spiral arm basin structure. 

Eigure 5 also reveals that destination times diverge on the basin boundaries (see 
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Figure 4. (A) Destination map section in the {d, '0)-plane with s = 0.56625 < Sc, for 

SD (red), DS (blue) and SSq (yellow) states. When s increases, the basin boundaries 
perform a spiraling motion as indicated by arrows. (B) The logarithmic times logT 
to destination reflect the structure of the destination map in (A). Times peak at the 
interface boundaries between SSq and SD/DS regions (see also Figure 5). Parameters 
are (A,/3) = (0.1,0.025). 



Figure 5. Times T to reach e-neighborhood of fixed points SSq/DS/SD; trajectories 
start from = (0.56625,0,7/;) for —tt < t/; < tt, A = 0.1,/3 = 0.025 (straight 

line on Figure 4). Average destination times T grow like a power law as t/; ^ Ttt, 
thus basin structure is self-similar around (d, t/;) = (0,7r). T diverges at the boundary 
between SD/DS and SSq basins (inset), since these trajectories lie on stable manifolds 
leading to saddle points on invariant manifolds Si,S 2 ^ 
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inset), which is explained by the fact that these boundaries form separatrix sheets: these 
are the two-dimensional stable manifolds emanating from the saddle chimeras on Si and 
S 2 , originating in the saddle-node bifurcation that gives birth to chimeras, see Figure 2 
and [41], Numerical continuation of those sheets (see Figs. 6 and E4, and Appendix 
E.7) displays the same twisting motion as seen in Figure 4. 


4. Analysis 

A complete analysis of the basins for the entire phase space is difficult to achieve, but 
the main features of the basin structure have their origin in the invariant rays Rq and 
about which a perturbation analysis can be made for small A and (5 (asymptotic results 
remain qualitatively in fair agreement for parameter values further off the origin). 


4 . 1 . Perturbation analysis around the invariant ray Rq 


We consider the coupling constants (/r, n) to be perturbed from global coupling {A = 0) 
by setting /r = |(1 -|- A), v = |(1 — A) as in [41]. We then make the perturbative 
approximation that xf,d,l3 and A are all small and of the same order (while keeping in 
mind that A > 2/5 is required for the existence of a chimera state in this limit). This 
means that near Rq, we make the ansatz [41]: 

if = + 

d = die + 0{e‘^), 

/3 = (die + 0{e‘^), 

A = Aie + Oie^) . 


After making a change of variables x = di, y = and then a second change of 

variables x = rcos{0) and y = rsin(6*) we find the following equations: 


^ = - sin(26») -F /5i^ sVe, 

^ ^ (1 - - 1 [1 + cos( 2 d)s 2 ] Aie, 

ds \ , 


/5ie. 


(5) 

( 6 ) 
(7) 


dt 2 

Note that the derivatives of r and s are both order e while the derivative of 9 is order 1 
(when s is not close to 1). Thus 9 evolves on a fast time scale while r and s evolve slowly. 
We may therefore use the method of averaging on the higher order terms involving 9 to 
simplify these equations to 

dr 2 

Tt = 

dO \ / 9 A \ 


II The change of variables is chosen so that that the spiraling cycles become circular in shape. 


( 8 ) 

( 9 ) 
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dt 2^^^ s)/3ie. (10) 

This reveals a couple of properties. First of all, as expected, solutions will be spirals 
around the Rq manifold. The radius r goes to 0 as t cx) (for the truncated equations). 
Thus, trajectories slowly converge toward the Rq manifold until the approximations 
break down when ^ = 0 and higher order terms become significant. In other words, 
the i^o-manifold is weakly attracting. The frequency of rotation is a; = |(1 — — Aie). 

So, as s —>■ Sc ~ \/l — A, the rotation frequency dO/dt —> 0. This is here referred 
to as the Sg-plane. In this plane, the trajectories cease to have a spiral character and 
instead begin to separate and evolve toward the fully synchronized state or the DS or 
SD chimeras. 

Note that there is an alternative way that the “critical plane” could be defined. 
Setting dO/dt = 0 in Eq. (6) yields a minimum s solution Sc = 1 — Aie + O(e^) 
(possible only for particular 0 values), thus Sc ~ 1 — ^ is the smallest value of s for 
which rotation about ray Rq may stop. The difference between the two expressions 
Sc ~ Vl ~ A ~ 1 — A/2 and Sg ~ 1 ~ ^ comes from whether averaging has been applied 
or not. In the former case (averaged equations), rotation about Rq stops on average 
over all 0; in the latter case, rotation about Rq stops only for some particular 6. 

By symmetry, if a trajectory originating at (s, d, ip) converges to the SD state, the 
trajectory originating at (s, —d, —ip) must converge to the DS state. Therefore, the 
position relative to a separating boundary in the Sc plane determines the final state. 
Numerical integration confirms that trajectories converging to SD and DS chimeras 
form opposing sides of a positively oriented double helix centered on Rq (red and blue 
in Figure 2 and Supplementary Video 2). 

The “rotation” of the basin boundary as s increases along the Rq manifold (indicated 
symbolically in Figure 4 A) can also be understood analytically in the perturbative limit 
close to Rq. Equation (10) can be solved explicitly to get 


sit) = So j ^s^ + (l-s^)e-<^At . 


( 11 ) 


Taking d = (1 — s^)/2 to lowest order (from (9)), we can substitute in for s(t) and then 
integrate from t = 0 to t = Rrit to approximate the total angle change about the Rq 
manifold over the course of the trajectory. Here tcrit represents the time at which the 
trajectory s{t) reaches the critical plane s = Sc- 


t, 


crit 


= T|n 

epi 


1-.2 ( 2 - A )^ 


A{4-A) 

Integrating 6 gives a total angle change 


= 

epi 


(1-Si 




(1 - 


2')„2 

cJSq 


A6 = —— In 
epi 


1 


So 


= /3-Mn(l-A/2)-/3-Mnso • 


-1 


( 12 ) 


(13) 


(This can also be written as Ad = ln(sc/so), and then this expression is valid for 
either definition of Sg.) Thus, the boundary angle is proportional to ln(l — A/2) — 
ln(s), yielding a rotation rate of {/3s)~^ as the section plane s varies uniformly. 








Basins of Attraction for Chimera States 


10 


Since the angle of a trajectory at the critical plane determines the basin the 
trajectory belongs to, the appearance of the basin boundary in a section plane orthogonal 
to the ray Rq is just a line with angle proportional to A6. 


4-2. Perturbation analysis around the invariant ray 

We can perform a similar analysis around the ray by a similar ansatz: 

Ip = TT + ipie + 0{e^), 
d = die + 0{e‘^), 

P = Pie + 0{e‘^), 

A = Aie + 0{e‘^) . 

This time we make an (analogous) change of variables x = di, y = pi and 

again convert to polar coordinates. This yields 

dr 
dt 


-Pi (1 cos (2d)) - 


1 — 


1 


sin (2d) 


re 


(14) 


dt 2 


+ 


I , '^o (l - COS (2d)) + l-pi sin (2d)s^ 

1 + s"' 2 


ds 

dt 


-s (l - s^) PiAi + 


1 + 
1 — 


sin(2d) 


(15) 


Again we find that the derivative of 9 is larger than the other derivatives, allowing us 
to reduce the system to 


dr 

dt 

^Pire, 



de 

-^Vl - + 


ll - s' 

dt 

2^ 

J 1 + s 

ds 

dt 

is(l-s")A7li 

eP 




(16) 

(17) 

(18) 

Similar to the results near Rq, this analysis reveals a spiraling motion around the 
Rt^ manifold, but here the radius diverges exponentially. In contrast to the previous 
case, three distinct time scales are present: the derivative of d is an order of magnitude 
larger than the derivative of r and two orders of magnitude larger than the derivative of 
s. This means that the rotation around the manifold and the radial divergence away 
from the manifold occur more quickly than the translation along the manifold. In other 
words, each trajectory (and consequently any basin boundary) winds around within 
a plane with approximately fixed s (see cross section in Figure 4 and Supplementary 
Video 3. Rotation around the manifold occurs at a faster rate than divergence away, 
and both occur faster than translation along the manifold. 


^ The change of variables is chosen so that that the spiraling cycles become circular in shape; the 
particular shape differs from the one near the i?o-nianifold. 
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We can understand two qualitative aspects of the basin boundaries seen in numerics: 
(1) the basin boundaries are linear in the Sc-plane near Rq, and (2) the basin boundaries 
have spiral shape near 

The invariant ray Rq is surrounded by the SSq basin, and as we move further 
away form Rq, we enter the SD and DS basins, respectively (Figs. 4 and E3, FlC). 
The boundaries (separatrices) between SSq/SD and SSq/DS basins are the two stable 
manifolds leading to the SD/DS chimera saddles. The relative width of the SSq and 
SD/DS basins varies with parameters A and (5 (Fig. E2), and the SSq basin is thin close 
to the origin, and for smaller (A,/5). 

To understand the basin shapes near Rq, it is helpful to consider trajectories 
generated by points along a straight line orthogonal to Rq. We will thereby use the 
asymptotic results (5)-(7) derived in the previous section valid for smaller A and {3. 
Consider the set of points along a line segment parameterized by k with |A:| < 1 where 
r = krQ, 9 = 9q and s = Sc- The trajectories generated by integrating the equations 
with initial conditions along that line segment will intersect the plane perpendicular 
to Rq at some later time at s = Sc + where 5 <C e. By symmetry, if a point along 
the line with k > evolves toward the DS chimera, then the corresponding point with 
/c < 0 will be mapped to the SD chimera. Similarly, if a point with k > 0 is mapped 
to the synchronized state, the corresponding point with A: < 0 will also be mapped 
to the synchronized state. Suppose Tq <C S, so that all points along the line segment 
are chosen to be arbitrarily close to the Rq manifold. According to Eqs. (5)-(7), d9 is 
independent of r, so the images of these points in the plane Sc + will remain collinear. 
Now dr = 0{S‘^) and ds = S, so not only will the points along this segment remain 
collinear; as s increases, they will also remain arbitrarily close to each other and to 
Rq. For at least one particular choice of 9q, this line segment will reach the surface of 
the cylinder in a direction tangent to the intersection of the invariant surfaces Si and 
S 2 . As discussed above, this intersection is itself an invariant manifold, and thus the 
entire line segment will be mapped to the synchronized state, being the only attractor 
on the manifold. Thus there exists a line segment in the Sc plane that lies in the basin 
of attraction for the SSq state and that separates the basins for SD and DS chimeras. 
This suggests that, at least sufficiently close to Rq, the basin boundary between SD and 
DS chimeras must be linear. 

Near R.,^ the picture is different. Because there are three distinct time scales, with 
the evolution of both 9 and r faster than the evolution of s, trajectories initially close 
to generate spirals within a fixed pfane perpendicufar to s — this is the origin of the 
spiraf shape of the basin structure near R^^. 

These quafitative arguments can be made rigorous in the fimit where the SSq basin 
becomes a set of measure zero (infinitesimal thickness). The basin boundaries 
(separatrices) are visualized in Fig. 6 (also see Fig. E3 for a close up of the basins 
near Rq). 
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Figure 6. Separatrix surfaces continued from the SD and DS saddle points on the Si 
(red) and S 2 (blue) manifolds, respectively, shown from different view angles (A, B, 
C, D). Continuation is performed as described in the text for A = 0.1 and j3 = 0.025. 
Black and white dots denote stable and unstable fixed points, respectively. 


6. Control strategies 

Determination of the structure of the basins of attraction for this system naturally 
invites the question of whether we can “control” the system. This can mean several 
things, among them: (1) can we intervene during an initial transient so as to direct 
the system to a desired equilibrium; (2) can we perturb the system to move it from 
one stable equilibrium to another; (3) can we stabilize an unstable equilibrium. Each 
of these questions can also be examined with the goal of finding an “optimal” strategy 
of some kind, where optimality is usually defined as minimizing some aspect of the 
intervention. To answer questions (1) and (2), knowledge of the basin structure in the 
thermodynamic limit N ^ 00 is clearly useful, at least for sufficiently large . 

A full exploration of control and intervention strategies is beyond the scope of this 
paper, but as a demonstration of the power of our approach, we have performed a simple 
experiment. We wish to take a system at equilibrium in the DS chimera state, and to 
perturb it sufficiently that it goes to a different equilibrium (either SD or SSq). We 
restrict ourselves to finite perturbations of the form 61 O], + Q, k = 1,, 

where Q quantifies a uniform phase shift to all oscillators in the synchronous group 
(population 2). Since the perturbation leaves the system state on the invariant DS 
manifold, the final state may change from DS to SSq. 
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In the thermodynamic limit, this is equivalent to holding pi and p 2 = 1 constant 
while perturbing xp via xp xp — Q. Thus we expect the minimal required perturbation 
Qmin to be determined by the size of the restricted DS basin of attraction (restricted to 
the surface p 2 = 1, the top surface of the cylinder shown in Figure 2). 



Figure 7. Minimal uniform perturbation of synchronous oscillator phases Qmin needed 
to escape from DS chimera state equilibrium. Dashed line indicates asymptotic value 
for —> oo. Dots indicate results from numerical experiments at each fixed value. 

Figure 7 shows the expected asymptotic value of Qmin determined from our study 
of the basins of attraction presented here, as well as the Qmin values determined via 
numerical experiment for finite values of . The precise threshold varies slightly 
depending on the initial phases but asymptotes to the value observed in the continuum 
limit. This good agreement confirms that (1) our knowledge of the basins of attraction 
in the thermodynamic limit can indeed inform control strategies for chimera states, and 
(2) insight into the finite system is indeed gained by analysis of the thermodynamic 
limit. 

Our analysis of the thermodynamic limit suggests that switching from the SSq state 
to a DS (or SD) chimera requires a perturbation to pi (or P 2 ), because pi = p 2 = 1 is 
an invariant manifold. Hence, while a uniform phase shift 1 —>■ + Q) will perturb 

xp^ it will not be sufficient to desynchronize one of the two populations. Instead, a 
nonuniform phase shift that decreases the value of pi (or P 2 ) can accomplish the desired 
switching behavior. 

Finally, we note that the control strategy presented here is quite naive. It is likely 
that more “optimal” strategies exist in the sense that the perturbation magnitudes could 
be reduced. 
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The probability that a random initial condition evolves to a chimera state, while 
important for real-world applications, has not been a frequent topic of investigation [17, 
25]. Here, we have provided a detailed mathematical analysis unveiling the basin 
structure for a very simple system with two populations, allowing for insight into the 
chimera’s relative rarity. It remains to be seen whether similar efforts applied to other 
models such as neuronal or pulse-coupled oscillators (where reduction methods have 
become available only very recently [50,51]) will bear fruit, and how basin structure in 
those systems will compare. 

Oscillators on a ring with finite-range coupling exhibit chimera states that are (very 
long) transients with chaotic dynamics [52]. However, extensive computational analysis 
of the finite system (1) displays no such transient behavior while chaotic behavior is 
absent [44,53]; this difference in dynamic behavior (along with others) seems to be 
related to differing coupling topology [44]. Moreover, very small oscillator systems with 
= 2,3,4 are shown to display asymptotic stability of chimera states [54]. 

Sampling the immense initial state space associated with the case of N‘^ < oo would 
be a burdensome task. Our analysis was facilitated by considering oo, allowing 

us to focus on the low-dimensional order parameter dynamics on the OA-manifold [47] . 
While the higher dimensional dynamics off this manifold poses a challenge in its own 
right [30] , the continuum theory allows us to gain useful insight by mapping the discrete 
to the continuous order parameter, ~ (identity for oo). Though 

bifurcation boundaries may blur for very low and the finely filigreed basin boundaries 
near R.,^ may break down, general basin structures will look similar even for moderate 
N^. 

The stable manifolds of the saddles near SSq divide phase space into simply 
connected basins of attraction (see Figs. 6, E4 and Supplementary Video 4). Basins 
near the fully synchronized (SSq) and chimera (SD/DS) states are simple in structure 
and relatively large"*", resulting in robustness to perturbations (see also Fig. 7). The 
twisting motion around the invariant rays i?o, 7 r, however, yields a complex basin structure 
which we explain analytically. As one approaches the basin density diverges, and 
basins become locally intermingled [55]: perturbations in that region affect the fate of 
a trajectory drastically. 

Continuum theory allows us to construct initial phase densities leading to chimera 
states via ^^{0, ^) = ^ [l + ("ggg Appendix A). If instead initial 

phases 0^ are sampled uniformly on [—7r,7r], the probability distribution for R„ is 
unimodal with mean ~ l/Vand variance ~ 1/N'^. The probability distribution 
for $ 0 - is invariant and uniform; this observation combined with the intermingled basin 
structure near R-j^ explains why in practice random initial phases can lead to both 
chimera and fully synchronized states. Thus, this implies in particular that chimera 

"*■ Local basin volumes of chimeras presumably scale like |xsd,ds “ XsaddleI, where x G denotes 
coordinates of stable (SD/DS) and unstable (SADDLE) chimera states in the reduced phase space. 


Basins of Attraction for Chimera States 


15 


states, given random initial conditions, are not always rare occurrences (depending on 
parameter values, see also Fig. 3). 

The results we have presented focus on the case of two populations. While two 
populations allow for multi-stability between the fully synchronous state (SSq) and 
two symmetrically equivalent chimera states (SD and DS), generalizations of such 
hierarchical structure to n > 2 populations [45,56] make accessible larger configuration 
spaces of size 2” by variation of the synchronization-desynchronization patterns. One 
may wonder if any of the structures studied here retains relevance in cases of more than 
two populations [46,56]. It can be shown that the invariant hyper-ray corresponding to 
Rq defined hj = p and = 0 exists for n > 2, and that the flow on this ray is p ^ 1. 
This suggests that the phase space is skeletonized similarly and a similar analysis may 
be feasible—a task left for a future study. 

Biologically, systems with multiple coexisting chimera attractors have been 
proposed to describe ‘metastable dynamics’ required to modulate neural activity 
patterns [57], or to encode memory. Indeed, localized dynamical states are directly 
related to function in neural networks [58,59]; localized synchrony has been widely 
studied in neural field models as bump states [51,60,61], and are phenomenologically 
similar to chimera states. It is worth noting that chimeras occur in models of neural 
activity [62-66]. 

Implementations of chimera states could likely be achieved in micro-(opto)-electro- 
mechanical oscillators [67, 68] where synchronization patterns may have technological 
applications. Conversely, as power grid network topologies evolve to incorporate growing 
sources of renewable power, the resulting decentralized, hierarchical networks [69, 70[ 
may be threatened by chimera states, which could lead to large scale partial blackouts 
and unexpected behavior. 

The potential for applications—or threats—makes the dynamic re-configuration 
and switching between chimera configurations (possibly modulating functional 
properties of the underlying oscillator network) particularly relevant [71]; applications 
that modulate functional properties can only be achieved using detailed knowledge of the 
basin structure. As a test of principle, we successfully implemented a simple algorithm 
to move the finite oscillator system between different equilibria, demonstrating that an 
understanding of the basins of attraction for the oo system has value. We hope 

that future work will explore the construction of efficient control strategies to stabilize 
or prevent chimera states, with applications across many fields. 
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We consider the Kuramoto-Sakaguchi model with non-local coupling between n 
populations [41,42] 

n N-' 

(A.l) 

cr'=l ^ 1=1 

where Of is the phase of the fcth oscillator k = belonging to population 

a = 1,... ,n. To facilitate comparison with previous work [41,42,45], we consider the 
case of symmetric coupling with The phase lag parameter a tunes between 

the regimes of pure sine-coupling {a = 0) and pure cosine-coupling (a = 7r/2). In what 
follows, we introduce the re-parameterized phase lag parameter /3 = 7r/2 — a, since for 
this type of system chimeras emerge in the limit of cosine-coupling [41,45], i.e. /5 —7- 0. 
To make further progress, we consider the thermodynamic limit, i.e., the case of 
oo oscillators per population. This allows for a description of the dynamics 
in terms of the mean-field order parameter [47-49]. Eqs. (A.l) then give rise to the 
continuity equation 

dP d 

^ + (A.2) 

where f'^{0,t) is the probability density of oscillators in population a, and v'^{0,t) is 
their velocity, given by 

n 

t) = U! 'y ^ Kaa' 

a=l 

Here we have dropped the superscripts to simplify notation: 9 means 9^, and 9' means 
9'^'. Following Ott and Antonsen [47,48], we consider probability densities along a 
manifold given by 

^ + + (A.4) 

where * denotes complex conjugation and acr(t) is given by 

aa{t) = j e^^r{9A)d9. (A.5) 

Defining a„{t) = p„(f)e''‘^'^d) where and represent mean-field order parameters, the 
governing equation can be reduced to the system [28,41,45] 

1-^2 

Pa = - J] -faS (A.6) 

( t '=1 

fa = >^ - - —^ J] K^^fp^j cos + (3). (A.7) 

a'=l 

The Ott/Antonsen manifold, in which the Fourier coefficients /n(t) of the probability 
density / satisfy fnif) = is globally attracting for a frequency distribution 


f^'{9'A)d9'. 


(A.3) 
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with non-zero width A [48]. For identical oscillators (A = 0), the dynamics for the 
problem (with n = 2 populations) can be described by reduced equations using the 
Watanabe/Strogatz ansatz [72], as shown in Pikovsky and Rosenblum [30]; the authors 
showed that Eqs. (1) may also be subject to more complicated dynamics than those 
described by the Ott/Antonsen ansatz. Studies by Laing [26, 73] investigated the 

dynamics using the Ott/Antonsen ansatz for n = 2 populations for the case of non¬ 

identical frequencies and found that the dynamics for sufficiently small A is qualitatively 
equivalent to the dynamics obtained for A = 0. It is therefore justified to discuss the 

dynamics for A —?• 0 representing the case of nearly identical oscillators using the 

Ott/Antonsen reduction. 


Appendix B. Governing equations for two populations 


We restrict our attention to the case of n = 2 populations. Accordingly, we define 
the coupling parameters Ku = K 22 = /i and K 12 = K 21 = n; by rescaling time 
we can eliminate one parameter so that 1 = /i + u without loss of generality The 
remaining parameter is redefined via A = /j. — n, expressing the disparity of coupling 
between the two neighboring populations. By virtue of the translational symmetry, 
(j)„ ^ -\- const., the dynamics of the system is effectively three dimensional. We 

introduce the angular phase difference ^ = 01 — 02 of the order parameter, and the 
resulting governing equations become 


Pi = 
P 2 = 



[HPi sin 0 -F z/p2 sin (0 - 0)], 



[/ip2 sin 0 upi sin (0 0)], 


(B.l) 

(B.2) 


0 


1 + P 2 
2p2 

2pi 


[PP 2 cos 0 upi cos (0 0)] 
[ppi cos 0 UP 2 cos (0 - 0)] . 


(B.3) 


where the state variables lie in the domain D = {(pi,p2,0) G IR^jO < Pi,P 2 A 1,0 G 

[-7r,7r]}. 

To investigate the basins of attraction, it proves useful to express the dynamics in 
terms of the sums and difference of the order parameters, i.e., we define 


5 — 2 ^Pi + P2), 


(B.4) 

(B.5) 


and 0 as above. These variables belong to the domain defined by 0 G [—vr, tt], s G [0,1] 
and d G [—a, a] with a{s) = | — || — s| (the back-transformation is pi = s -|- d and 
P 2 = s — d, without the factor of 2.). The governing equations are then expressed as 




3d^ — s‘^) + ^{l + d'^ — cos 0] sin0 
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-|- vd{l — d^ s^) cos j3 sin'll. 

(B.6) 

-{(i[/r(l — d^ — 3s^) — 1^(1 — -|- cos 'ip] sin (3 


— ps{l + d^ — s^) cos j3 sin-^}. 

(B.7) 

(d^ — — 2 ds cos /3[/i(d^ — s^) -\- p cos ip\ 


-[- [s^ -|- d^ -|- -|- d^(l — 2 s^)\p sin (3 sixiip}. 

(B.8) 


Eqs. (B.1)-(B.3) or Eqs. (B.6)-(B.8), respectively, are invariant under the 
transformation E : {pi,p 2 ,'ip) i-)- {p 2 , pi, —f’), corresponding to interchanging the two 
oscillator populations. More generally, the change of parameters Yi \ (3 ^ {5 + tt reverses 
time in the governing equations, thus inverting flow and stability properties in phase 
space. This is also valid for the more general case of n > 2 equations, i.e., for Eqs. (A.6) 
and (A.7). 

Appendix C. Invariant manifolds (IMs). 

Two-dimensional invariant manifolds. Letting pi —>■ 1 in Eqs. (B.1)-(B.3) leaves pi 
invariant, i.e. pi = 0. The same holds true for p 2 by symmetry. Thus we find two two- 
dimensional invariant surfaces, corresponding to the top and side surface of D, defined 
by 

^cr = {{Pl, P2, G A I Pct = 1}) 

where cr = 1,2 refers to the SD, DS manifolds, respectively. The dynamics in one 
manifold is identical to the other via the symmetry operation defined by operator E, 
see main text. The dynamics on these IMs is analyzed in [41]. 

One-dimensional invariant manifolds. Our numerical investigations indicate the 
presence of an invariant manifold at V' = 0,7r with pi = P 2 - Substituting these values 
into Eqs. (B.1)-(B.3), we get 

s = - sin/3 • s(l — s^), (C.2) 

d = ^ = 0. (C.3) 

The first equation implies that any initial point on the rays with d = 0 and = 0,7r 
remains there for all times; if 0 < /3 < tt, the trajectory moves towards the SSq attractor 
according to (C.2). Thus, two invariant rays exist, defined via 

R 4 , = {(Pi, P2, e D\pi = p 2 and = f} (C.4) 

with 0 = 0, TT. 

Note that another one-dimensional invariant manifold S'12 is defined as the 
intersection S\ 0 S 2 , and any initial point with s = 1 on S'12 will therefore always 
end up in the SSq state. 
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Appendix D. Fixed points. 

Fixed points on Si^2- The fixed points in the S' 1,2 manifolds are the SD, DS chimera 
states and fully synchronized SSq states that are discussed in detail in [41]; note that 
since 5 'i 2 is an invariant manifold, there must be another fixed point in addition to SSq 
contained in it, with opposite stability: this source is found at {pi,p 2 ,'ip) = ( 1 , 1 ,tt), 
which we refer to as SS-„-- Figure 2 illustrates how trajectories nearby are repelled from 
the ray On S' 1 ^ 2 , stable chimera states are born through a saddle node bifurcation, 
and undergo a Hopf bifurcation for sufficiently large disparity values A so that po- < 1 is 
oscillatory; the associated limit cycle is destroyed in a homoclinic bifurcation with even 
larger A. 

Chimera states. In addition to the in-phase (pi = p2 = 1 and = 0) and anti-phase 
(pi = P 2 = 1 and = w) equilibrium points, there are also three equilibrium points 
with p 2 = 1 and pi 7 ^ 1 (and three analogous fixed points with pi = 1 and p 2 7 ^ 1 ) [28]. 
These equilibrium points represent chimera states. Numerics suggest that two of these 
equilibrium points occur near = 0 and one occurs near = tt. Using an ansatz 
motivated by these numerical results, we find that these equilibrium points satisfy the 
following scaling relationships (where A = p — u and p + p = 1: 

i) Stable Chimera near = b) (DS): 


P ~ APi, 



a) Unstable Saddle Chimera near ip = 0 (UC): 
P ~ APi, 



Hi) Unstable Chimera near ip = tt (UC): 
P ~ A^Pi, 
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These relationships are useful when trying to solve for the precise fixed point 
locations numerically and for approximating their stable and unstable manifolds in order 
to deduce the basin boundaries. We note that chimera states asymptotically approach 
either SSq or SS^ as {A,P) (0,0). 

Origin. The origin is an unstable fixed point, as can be seen by linearizing for small 
Pi and p 2 in Eqs. (B.1)-(B.3). 


Other fixed points. Here we ask whether there are any fixed points off the invariant 
manifolds 5'i,2: i.e., are there any fixed points with pi,p 2 ^ *S'i ,2 where no population 
is completely synchronized? Together with Eqs. (B.1)-(B.3) ,0< Pi, P 2 < 1 implies the 
following conditions: 

0 =/rpi sin/? + z^P 2 sin (/? —'^), (D.l) 

0 =/ip 2 sin/? + sin (/? +//;), (D.2) 

0 = pi(l + pI) [p,p 2 cos/3 + upi cos (/3 + f;)] 

- P2{1 + Pi) [ppi cos/3 + z/p 2 Cos (/3 - f;)]. (D.3) 

We know that /3 —?• /3 + vr reverses time, so we can w.l.o.g. restrict our attention to 
0 < /3 < TT. When /3 = 0, tt, the first two equations are satisfied if ?// = 0, tt. The third 

equation yields the solutions p 2 = ±pi and p 2 =-where 

only the first branch lies in 0 < pi ,2 < 1- 

For all other cases, let us consider the equations by introducing K = p/u > 1 and 

Prel = P2/Pi- 

0 = [cos {'tp)Prei + K] sin/3 — p„i cos /3 sin fj, (D.4) 

0 = [Kp^^i + cosfj] sin/3 + cos/3sin//;, (D.5) 

0 = - [2pLp? + (pIi - l)Pi] sin 1 /; sin/3 (D.6) 

+ [(1 - pDpi cos//^ + (p2, - l)Kp,,^pl] cos/3. 

We note now that p^^i > 0 and sin /3 > 0 by assumption and we can eliminate these 
expressions as follows 

cos /3 sin 1/3 cos (//^) p^ei + K 


sin /3 

cos /3 sin Ip 
sin /3 


Prel 

= — Kp^^i — COSp!. 


(D.7) 

(D.8) 


Equating (D.7) and (D.8), it follows that a fixed point with 0 < pi, p 2 < 1 and 0 < /3 < tt 
can only exist if 

0 = TCp^gj + 2 cos {tp) p„i + K, (D.9) 

which has the solutions 

cos Ip ± cos^ (ip) — K'^ 


Prel 


K 


(D.IO) 
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However, by assumption, we have K > 1 and the solutions are complex. Therefore, 
even if we find real values 'tp,pi that satisfy the third fixed point equation (D.6), there 
will be no real solutions for p,^i, and thus also not for p 2 - Therefore fixed points in the 
interior of the domain can be excluded when (5 ^ mr where n is an integer. 

Appendix E. Numerical Analysis 

Appendix E.l. Probabilistic measure of basins of attraction 

In order to obtain an estimate of the sizes of the basins of attraction of the equilibria 
of (B.1)-(B.3), we selected 1000 random initial points {pi,p 2 ,'ip). Eqs. (B.1)-(B.3) were 
then integrated for a sweep of parameter values of 0.01 < A < 0.49 with increments of 
0.01 and 0.005 < P < 0.245 with increments of 0.005 until a final state was detected. 
The contour plot in Figure 3 displays the fraction of those trajectories with final states 
near a chimera state. 

It should be noted that the numerical experiment above assumes that pi, p 2 , and 
are uniformly distributed. For systems with a finite number of oscillators N'^ in each 
population, the expected value of the order parameter value p^- is 0{1/'/N^). Hence, 
the probabilities computed using the above scheme should not be interpreted as the 
probability that a state with randomly selected initial phases 9^^^ would evolve toward 
a chimera state. Instead, they represent the size of the basins of attraction of the chimera 
states relative to the size of the basin of attraction of the fully synchronized state in the 
continuum limit N‘^ oo. 

Appendix E.2. Destination Maps 

Simulations for a given initial condition were carried out until trajectories to a fully 
synchronized (limit point, LP) or a stable (LP) or breathing chimera (limit cycle 
LC) occurred. The detection of these three types of states was carried out in two 
steps, described below. Integration of Eqs. (B.1)-(B.3) or (B.6)-(B.8) were carried 
out in Matlab^^ using the ode45 solver routine with event detection (see below) 
on a high performance computation cluster, with a relative error tolerance of 10“®. 
Algorithms below are outlined for (pi,p 2 ,l/’)-coordinates; analogous detections for 
(s, d, 'i/^j-coordinates are carried out by applying the related coordinate transformations. 
Below, dpi ,2 ~ Pi, 2 dt and dp ^ dt denote the approximate differential values evaluated 
by the o.d.e. integrator at discrete time steps. 

Appendix E.3. Simple convergence test (Event A) 

This simple test was used to detect the type of state is asymptotically achieved. 
Integration was stopped by an event detection algorithm solving for roots of: 

• LP detection: v = [dp1 + dp 2 + — S: convergence to any LP (in any direction). 
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• Convergence to LC on v = [(pi — 1)^ + dp^] — S, passage through p 2 = 0 (= 1 
cycle), positive direction. 

• Convergence to LC on S 2 '. v = [{p 2 — 1)^ + dpf\ — 5, passage through pi = 0 (= 1 
cycle), positive direction. 

A convergence tolerance of 5 ~ 10“® was chosen. 

Appendix E.f- Estimating time to attractor 

The following algorithm is adopted for obtaining estimates for the time to reach the 
attractor, T, i.e., the traveling time from initial to end condition. When these times are 
not of interest, the previous scheme is preferred due to significant gains in computational 
speed. 

(i) Integration is carried out until v = [dpi + dpi + — 5 crosses a zero (Event 

B, LP detection). 

(ii) If the integration fails to detect a fixed point, the algorithm enters a loop of max. 
100 iterations, where: 

i. ) Integration is carried out to detect k = 1,..., 10 events of type Event A, 
limit point and limit cycles). Periods of limit cycles and event states {pi,P 2 , 'f’)\t=tk 
are stored. 

ii. ) Test for convergence to limit point or limit cycle: \\{pi, P2,'4’)\t=tk ~ 

(pi,P2,^)|t=tfc+ill < Cc with Cc ~ 10“^ 

hi.) Exit loop when a LP or LC is detected or 100 iterations are carried out. 

(hi) If LC or LP is detected, the final state is detected as explained above. Otherwise, 
failed convergence is stored as a failed end state. 


Appendix E.5. Destination maps in the Sc-plane 

Destination maps were calculated for (5 = 0.01,..., 0.125 at constant A = 0.2 
(Figure E2A) and for A = 0.08,..., 0.41 at constant /3 = 0.05 (Figure E2B). Saddle-node 
(SN) and homoclinic (HC) transitions are indicated. 

Appendix E.6. Destination map for small s 

In Figure E3, we display a sample destination map computed for small s for the purpose 
of demonstrating what the basin structure would look like if initial phases were chosen 
from a uniform random distribution. If, for example, = 50, the expected initial 
value of s would be close to 0.1. 

Appendix E.7. Numerical continuation of the basin boundaries (separatrices) 

The stable manifold of the saddle chimera defines the boundary of the basin of attraction 
of the corresponding stable chimera. By approximating this manifold, we can visualize 
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Figure El. Times to attractor. Parameter values for (A, B): A = 0.1,= 0.025 
at pi = p 2 = 0.56625 and for (C, D): A = 0.1,/d = 0.05 at pi = p 2 = 0.5. Final 
destinations are color coded in red for SD, blue for DS and yellow for SSq states. 


a 


d 



Figure E2. Destination maps for parameter sweeps in A (vertical) and (3 (horizontal), 
respectively. The maps are shown in the (d, '0)-plane at s = 5c = 1 — A. Parameters 
where saddle node and homoclinic bifurcations occur are denoted by SN and HC, 
respectively. 
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if 

Figure E3. Destination map in the {d,ip) plane with s = 0.1. Red indicates SD 
chimera, blue indicates DS chimera, yellow indicates SSq state. Parameters are A = 0.1 
and p = 0.025. 


which regions of the state space will evolve toward this chimera state, as shown in 
Figure 6 and Figure E4. The manifold can be approximated as follows: 

Step 1: Compute the two stable eigenvectors of the saddle chimera to obtain a 
local approximation to the stable manifold near the saddle chimera. (There are two 
stable eigenvectors vi and V 2 and one unstable eigenvector V 3 for the saddle chimera 
located at p. The stable eigenvectors define a plane tangent to the stable manifold.) 

Step 2 : Obtain a family of starting points xo(0) for continuation by making small 
perturbations off of the saddle chimera in every direction within the stable manifold. 
(Define a vector of angles 9 and a magnitude e. The family of starting points are defined 
by Xo(d) = p + e [cos(d)vi + sin(d)v 2 ].) 

In Figure 6 , we used 23 angles 9 between 0 and vr with the vectors Vi and V 2 chosen 
so that all of these perturbations led to relevant parameter values and a perturbation 
magnitude of e = 10“®. In Figure E4, 94 trajectories were used. 

Step 3: Integrate backward in time from each point until the trajectories reach 
Xi with a predetermined distance ||xi — Xo|| from the start point, and plot the surface 
containing these trajectories. (We used ode45 to integrate the equations and then 
interpolated to determine when the trajectories had reached the desired length.) 

In Figure 6 , the predetermined distance was set at 0.01 to obtain a high resolution 
near the manifolds. In Figure E4, the predetermined distance was set between 1 and 20 
(and no additional refinement was performed) in order to reduce data points for quick 
rendering, and in particular to enhance the visibility of the separatrices while reducing 
the total number of points displayed. This way the point families are equidistant in 
space, rendering an accurate picture of the separatrix surface in all regions. 

Step 4: The endpoints of the trajectories define a curve. Use evenly spaced points 
along the curve as new starting points and return to step 3 until enough of the stable 
manifold has been computed. 

In Figure 6 , we used a spacing of 0.01 near the saddle chimera and 0.05 once the 
trajectories had reached a distance of 0.2 from the saddle chimera. 
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Figure E4. Visualization of separatrix surfaces and trajectories. Points along the 
separatrix corresponding the DS chimera state are colored blue and points along the 
separatrix corresponding the SD chimera state are colored red. A: Continuation of the 
separatrix for the DS chimera state. Stable manifold and corresponding eigenvectors 
of SADDLE are shown solid (magenta) and dashed (green), and unstable eigenvector 
in yellow. Red dots indicate initial points from where the stable manifold (blue) 
was continued. B: Superposition of the two continued separatrix surfaces. The 
continuation in A and B is performed as described in the Appendix. C: Trajectories 
along the separatrix surfaces originating from SD and DS saddles points on the Si 
and S 2 manifolds, respectively. D: Trajectories (dotted) along the separatrix surfaces, 
continued from saddle chimeras. Parameters are (3 = 0.025 and A = 0.1 (A-C) or 
A = 0.2 (D). 


While this method yields satisfactory results for the problem at hand, we 
mention that more advanced and accurate continuation methods are available for the 
computation of the manifold, for an overview of such methods see [74]. 

Appendix F. Alternative coordinate representation 

In the main text, we chose to use the parametrization with pi,p 2 and '0, because 
this allows for visualization in a familiar cylindrical coordinate system, and because 
these coordinates have natural interpretations in terms of the distributions of phases 
in the finite oscillator system: pi and p 2 indicate the degree of synchrony in each 
population, and defines the mean phase difference between the populations. However, 
an alternative coordinate representation is possible that better reflects the symmetries 
inherent to the system, as is discussed here. 
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The equations describing the thermodynamic limit (B.1)-(B.3), before being 
transformed into polar coordinates, can be rewritten in terms of two complex amplitudes, 
Zk = PkC'^’^, k = 1,2, taking the form 

and the corresponding eqnation for Z 2 with interchanged indices. This system exhibits 
a rotational symmetry according to 

(zi,Z2) (zie*‘^,Z2e*‘^)- 

This symmetry motivates a reduced coordinate system 

7 = ^1^2 E C (F.2) 

h = |;2 i| 2 _ e M (F.3) 

with 0 < I7I < 1, and h E [—1,1], from which we recover the original variables with the 
intnitive meaning 

P2 = N|" = ^(\/<(^ + 4|7P-<() 

f; = arg (zi - Z 2 ) = arg c, 

provided that \zi\ and \z 2 \ are non-zero. This system is singnlar only a,t zi = Z 2 = 0 
and its geometry can be presented so that the symmetry of exchanging Zi and Z 2 is 
maintained, i.e., the reflection symmetry 

7^7. 

For this parameterization, the fully synchronized states SSo and SSjr are located at 
(7,h) = (1,0) and (7,(5) = (—1,0), respectively. The invariant rays Rq and are 
located on the same straight line given by (5 = 0 with Imy = 0. The invariant manifolds 
Si and S 2 are the two paraboloids defined via ±(5 = 1 — jyp. States of interest are 
then in the region enclosed by these two paraboloids. Sample trajectories are shown in 
Fig. FI. 
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Figure FI. Trajectories with initial conditions near Rq (A, close up in C) and 
near manifold (B), leading to the states SD (red), DS (blue) and SSq (yellow). 
Parameters are A = 0.1 and (3 = 0.025. The invariant surface manifolds Si and S 2 are 
colored red and blue, respectively. 
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